#----Basics----

rm(list = ls())

# Set working directory
setwd("~/Dropbox/Whitten_Kagalwala/TS/JOP/Acceptance/KW_Replication/UnitRoot_Independent")

# Load required libraries
library(tidyverse)
library(haven)
library(stringr)

# Load data
data <- read_dta('unitroot_independent.dta')
data$time <- 1:7

#----Reshape Data----
phi1_bias <- data %>%
  dplyr::select(c(time, starts_with("phi_bias"), starts_with("phi1_bias"))) %>%
  pivot_longer(phi_bias_pa:phi1_bias_adl33, names_to = "model", values_to = "phi1_bias") %>%
  mutate(model = gsub("^ph.*_bias_(.*)", "\\1", model))
phi1_type1 <- data %>%
  dplyr::select(c(time,matches("^phi_(.*)_t1$"), matches("^phi1_(.*)_t1$"))) %>%
  pivot_longer(phi_pa_t1:phi1_adl33_t1, names_to = "model", values_to = "phi1_type1") %>%
  mutate(model = gsub("^ph.*_(.*)_t1", "\\1", model))
phi1_power <- data %>%
  dplyr::select(c(time,matches("^phi_(.*)_power$"), matches("^phi1_(.*)_power$"))) %>%
  pivot_longer(phi_pa_power:phi1_adl33_power, names_to = "model", values_to = "phi1_power") %>%
  mutate(model = gsub("^ph.*_(.*)_power", "\\1", model))
phi2_bias <- data %>%
  dplyr::select(c(time, starts_with("phi2_bias"))) %>%
  pivot_longer(phi2_bias_adl22:phi2_bias_adl33, names_to = "model", values_to = "phi2_bias") %>%
  mutate(model = gsub("^ph.*_bias_(.*)", "\\1", model))
phi2_type1 <- data %>%
  dplyr::select(c(time,matches("^phi2_(.*)_t1$"))) %>%
  pivot_longer(phi2_adl22_t1:phi2_adl33_t1, names_to = "model", values_to = "phi2_type1") %>%
  mutate(model = gsub("^ph.*_(.*)_t1", "\\1", model))
phi3_bias <- data %>%
  dplyr::select(c(time, starts_with("phi3_bias"))) %>%
  pivot_longer(phi3_bias_adl33, names_to = "model", values_to = "phi3_bias") %>%
  mutate(model = gsub("^ph.*_bias_(.*)", "\\1", model))
phi3_type1 <- data %>%
  dplyr::select(c(time,matches("^phi3_(.*)_t1$"))) %>%
  pivot_longer(phi3_adl33_t1, names_to = "model", values_to = "phi3_type1") %>%
  mutate(model = gsub("^ph.*_(.*)_t1", "\\1", model))
b1_bias <- data %>%
  dplyr::select(c(time, starts_with("b1_bias"))) %>%
  pivot_longer(b1_bias_stat:b1_bias_adl33, names_to = "model", values_to = "b1_bias") %>%
  mutate(model = gsub("^b1_bias_(.*)", "\\1", model))
b1_type1 <- data %>%
  dplyr::select(c(time,matches("^b1_(.*)_t1$"))) %>%
  pivot_longer(b1_stat_t1:b1_adl33_t1, names_to = "model", values_to = "b1_type1") %>%
  mutate(model = gsub("^b1_(.*)_t1", "\\1", model))
b2_bias <- data %>%
  dplyr::select(c(time, starts_with("b2_bias"))) %>%
  pivot_longer(b2_bias_ds:b2_bias_adl33, names_to = "model", values_to = "b2_bias") %>%
  mutate(model = gsub("^b2_bias_(.*)", "\\1", model))
b2_type1 <- data %>%
  dplyr::select(c(time,matches("^b2_(.*)_t1$"))) %>%
  pivot_longer(b2_ds_t1:b2_adl33_t1, names_to = "model", values_to = "b2_type1") %>%
  mutate(model = gsub("^b2_(.*)_t1", "\\1", model))
b3_bias <- data %>%
  dplyr::select(c(time, starts_with("b3_bias"))) %>%
  pivot_longer(b3_bias_adl22:b3_bias_adl33, names_to = "model", values_to = "b3_bias") %>%
  mutate(model = gsub("^b3_bias_(.*)", "\\1", model))
b3_type1 <- data %>%
  dplyr::select(c(time,matches("^b3_(.*)_t1$"))) %>%
  pivot_longer(b3_adl22_t1:b3_adl33_t1, names_to = "model", values_to = "b3_type1") %>%
  mutate(model = gsub("^b3_(.*)_t1", "\\1", model))
b4_bias <- data %>%
  dplyr::select(c(time, starts_with("b4_bias"))) %>%
  pivot_longer(b4_bias_adl33, names_to = "model", values_to = "b4_bias") %>%
  mutate(model = gsub("^b4_bias_(.*)", "\\1", model))
b4_type1 <- data %>%
  dplyr::select(c(time,matches("^b4_(.*)_t1$"))) %>%
  pivot_longer(b4_adl33_t1, names_to = "model", values_to = "b4_type1") %>%
  mutate(model = gsub("^b4_(.*)_t1", "\\1", model))
lrm_bias <- data %>%
  dplyr::select(c(time, starts_with("lrm_bias"))) %>%
  pivot_longer(lrm_bias_pa:lrm_bias_adl33, names_to = "model", values_to = "lrm_bias") %>%
  mutate(model = gsub("^lrm_bias_(.*)", "\\1", model))
lrm_type1 <- data %>%
  dplyr::select(c(time,matches("^lrm_(.*)_t1$"))) %>%
  pivot_longer(lrm_pa_t1:lrm_adl33_t1, names_to = "model", values_to = "lrm_type1") %>%
  mutate(model = gsub("^lrm_(.*)_t1", "\\1", model))
autocorr <- data %>%
  dplyr::select(c(time,matches("^(.*)_autocorr$"))) %>%
  pivot_longer(stat_autocorr, names_to = "model", values_to = "autocorr_type1") %>%
  mutate(model = gsub("^(.*)_autocorr", "\\1", model))
ftest_y_t1 <- data %>%
  dplyr::select(c(time,matches("^ftest_y_(.*)_t1$"))) %>%
  pivot_longer(ftest_y_adl33_t1:ftest_y_adl22_t1, names_to = "model", values_to = "ftest_y_t1") %>%
  mutate(model = gsub("^ftest.*_(.*)_t1", "\\1", model))
ftest_x_t1 <- data %>%
  dplyr::select(c(time,matches("^ftest_x_(.*)_t1$"))) %>%
  pivot_longer(ftest_x_adl33_t1:ftest_x_adl22_t1, names_to = "model", values_to = "ftest_x_t1") %>%
  mutate(model = gsub("^ftest.*_(.*)_t1", "\\1", model))
ftest_x_t12 <- data %>%
  dplyr::select(c(time,matches("^ftest_x_(.*)_j1_t1$"))) %>%
  pivot_longer(ftest_x_adl33_j1_t1:ftest_x_adl22_j1_t1, names_to = "model", values_to = "ftest_x_t1") %>%
  mutate(model = gsub("^ftest.*_(.*)_j1_t1", "\\1", model))
ftest_y_power <- data %>%
  dplyr::select(c(time,matches("^ftest_y_(.*)_power$"))) %>%
  pivot_longer(ftest_y_adl33_power:ftest_y_adl22_power, names_to = "model", values_to = "ftest_y_power") %>%
  mutate(model = gsub("^ftest.*_(.*)_power", "\\1", model))
ftest_y_power2 <- data %>%
  dplyr::select(c(time,matches("^ftest_y_(.*)_j1_power$"))) %>%
  pivot_longer(ftest_y_adl33_j1_power:ftest_y_adl22_j1_power, names_to = "model", values_to = "ftest_y_power") %>%
  mutate(model = gsub("^ftest.*_(.*)_j1_power", "\\1", model))

#---- Manuscript Figures 1-7 Below ----

#----Figure 1, Granger & Newbold (1974) Replication----
b1_type1 %>% 
  filter(model == "stat") %>%
  mutate(gn = ifelse(time == 2, "1", "0")) %>%
  ggplot(aes(shape = gn)) +
  geom_hline(yintercept = 0.05, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = b1_type1)) +
  theme_bw() +
  labs(y = "Type 1 Error Rate", x= "T", title = "Granger and Newbold (1974) Replication") + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10)) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("1", "0"),
                     values = c(15, 0)) + 
  theme(legend.position = "none")
ggsave("GN1974.pdf" , width = 6, height = 4, units = "in")
ggsave("GN1974.eps", device = "eps", width = 6, height = 4, units = "in")

#----Figure 2, Residual autocorrelation in the static model----
data %>%
  ggplot() +
  geom_point(aes(x = time, y = corr_coeff), shape = 0) +
  theme_bw() +
  labs(y = expression("Correlation Coefficient"), x= "T", title = "Residual Autocorrelation in the Static Model") + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10)) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  theme(legend.position = "none")
ggsave("residauto_static_corr_coeff.pdf" , width = 6, height = 4, units = "in")
ggsave("residauto_static_corr_coeff.eps", device = "eps", width = 6, height = 4, units = "in")

#----Figure 3a, ADL(1,1) model type 1 error rates - x ----
b1_type1 %>% 
  filter(model == "adl") %>%
  ggplot() +
  geom_hline(yintercept = 0.05, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = b1_type1), shape = 4) +
  theme_bw() +
  labs(y = "Type 1 Error Rate", x= "T", title = expression("Type 1 error rate in the coefficient of x"[t])) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10)) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  theme(legend.position = "none") +
  ylim(0, 0.6)
ggsave("UnitRoot_Independent_b1_type1_ADL11.pdf" , width = 6, height = 4, units = "in")
ggsave("UnitRoot_Independent_b1_type1_ADL11.eps", device = "eps", width = 6, height = 4, units = "in")

#----Figure 3b, ADL(1,1) model type 1 error rates - Lx ----
b2_type1 %>% 
  filter(model == "adl") %>%
  ggplot() +
  geom_hline(yintercept = 0.05, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = b2_type1), shape = 4) +
  theme_bw() +
  labs(y = "Type 1 Error Rate", x= "T", title = expression("Type 1 error rate in the coefficient of x"[t-1])) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10)) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  theme(legend.position = "none") +
  ylim(0, 0.6)
ggsave("UnitRoot_Independent_b2_type1_ADL11.pdf" , width = 6, height = 4, units = "in")
ggsave("UnitRoot_Independent_b2_type1_ADL11.eps", device = "eps", width = 6, height = 4, units = "in")

#----Figure 3c, ADL(1,1) model type 1 error rates - Ly ----
phi1_type1 %>% 
  filter(model == "adl") %>%
  ggplot() +
  geom_hline(yintercept = 0.05, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = phi1_type1), shape = 4) +
  theme_bw() +
  labs(y = "Type 1 Error Rate", x= "T", title = expression("Type 1 error rate in the coefficient of y"[t-1])) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10)) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  theme(legend.position = "none") +
  ylim(0, 0.6)
ggsave("UnitRoot_Independent_phi1_type1_ADL11.pdf" , width = 6, height = 4, units = "in")
ggsave("UnitRoot_Independent_phi1_type1_ADL11.eps", device = "eps", width = 6, height = 4, units = "in")

#----Figure 4a, x Bias----
b1_bias$model <- factor(b1_bias$model, levels = c("stat", "diff", "pa", "gecm", "adl", "adl22", "adl33"), labels =  c("stat", "diff", "pa", "gecm", "adl", "adl22", "adl33"))
ggplot(b1_bias, aes(shape = model)) +
  geom_hline(yintercept = 0, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = b1_bias), position = position_dodge(width = 0.75)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 2, byrow = TRUE)) +
  labs(y = "Bias", x= "T", shape = "", title = expression("Bias in the coefficient of x"[t])) + ylim(-0.5, 0.5) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("stat", "diff", "pa", "gecm", "adl", "adl22", "adl33"),
                     values = c(0, 1, 2, 3, 4, 5, 6),
                     labels = c("Static", "Differences", "PA", "GECM", "ADL(1,1)", "ADL(2,2)", "ADL(3,3)")) +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  ylim(-0.1, 0.1)
ggsave("UnitRoot_Independent_b1_bias.pdf", width = 6, height = 4, units = "in")
ggsave("UnitRoot_Independent_b1_bias.eps", device = "eps", width = 6, height = 4, units = "in")

#----Figure 4b, x Type 1----
b1_type1$model <- factor(b1_type1$model, levels = c("stat", "diff", "pa", "gecm", "adl", "adl22", "adl33"), labels =  c("stat", "diff", "pa", "gecm", "adl", "adl22", "adl33"))
ggplot(b1_type1, aes(shape = model)) +
  geom_hline(yintercept = 0.05, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = b1_type1), position = position_dodge(width = 0.75)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 2, byrow = TRUE)) +
  labs(y = "Type 1 Error Rate", x= "T", shape = "", title = expression("Type 1 error rate in the coefficient of x"[t])) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("stat", "diff", "pa", "gecm", "adl", "adl22", "adl33"),
                     values = c(0, 1, 2, 3, 4, 5, 6),
                     labels = c("Static", "Differences", "PA", "GECM", "ADL(1,1)", "ADL(2,2)", "ADL(3,3)")) +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1)
ggsave("UnitRoot_Independent_b1_type1.pdf", width = 6, height = 4, units = "in")
ggsave("UnitRoot_Independent_b1_type1.eps", device = "eps", width = 6, height = 4, units = "in")

#----Figure 5a, Lx Bias----
b2_bias$model <- factor(b2_bias$model, levels = c("ds", "gecm", "adl", "adl22", "adl33"), labels =  c("ds", "gecm", "adl", "adl22", "adl33"))
ggplot(b2_bias, aes(shape = model)) +
  geom_hline(yintercept = 0, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = b2_bias), position = position_dodge(width = 0.75)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 2, byrow = TRUE)) +
  labs(y = "Bias", x= "T", shape = "", title = expression("Bias in the coefficient of x"[t-1])) + ylim(-0.5, 0.5) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("ds", "gecm", "adl", "adl22", "adl33"),
                     values = c(8, 3, 4, 5, 6),
                     labels = c("Dead Start", "GECM", "ADL(1,1)", "ADL(2,2)", "ADL(3,3)")) +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  ylim(-0.1, 0.1)
ggsave("UnitRoot_Independent_b2_bias.pdf", width = 6, height = 4, units = "in")
ggsave("UnitRoot_Independent_b2_bias.eps", device = "eps", width = 6, height = 4, units = "in")

#----Figure 5b, Lx Type 1----
b2_type1$model <- factor(b2_type1$model, levels = c("ds", "gecm", "adl", "adl22", "adl33"), labels =  c("ds", "gecm", "adl", "adl22", "adl33"))
ggplot(b2_type1, aes(shape = model)) +
  geom_hline(yintercept = 0.05, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = b2_type1), position = position_dodge(width = 0.75)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 2, byrow = TRUE)) +
  labs(y = "Type 1 Error Rate", x= "T", shape = "", title = expression("Type 1 error rate in the coefficient of x"[t-1])) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("ds", "gecm", "adl", "adl22", "adl33"),
                     values = c(8, 3, 4, 5, 6),
                     labels = c("Dead Start", "GECM", "ADL(1,1)", "ADL(2,2)", "ADL(3,3)")) +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1)
ggsave("UnitRoot_Independent_b2_type1.pdf", width = 6, height = 4, units = "in")
ggsave("UnitRoot_Independent_b2_type1.eps", device = "eps", width = 6, height = 4, units = "in")


#----Figure 6a, Ly Bias----
phi1_bias$model <- factor(phi1_bias$model, levels = c("pa", "ds", "gecm", "adl", "adl22", "adl33"), labels = c("pa", "ds", "gecm", "adl", "adl22", "adl33"))
ggplot(phi1_bias, aes(shape = model)) +
  geom_hline(yintercept = 0, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = phi1_bias), position = position_dodge(width = 0.75)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 2, byrow = TRUE)) +
  labs(y = "Bias", x= "T", shape = "", title = expression("Bias in the coefficient of y"[t-1])) + ylim(-0.5, 0.5) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("pa", "ds", "gecm", "adl", "adl22", "adl33"),
                     values = c(2, 8, 3, 4, 5, 6),
                     labels = c("PA", 'Dead Start', "GECM", "ADL(1,1)", "ADL(2,2)", "ADL(3,3)")) +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  ylim(-0.3, 0.1)
ggsave("UnitRoot_Independent_phi1_bias.pdf", width = 6, height = 4, units = "in")
ggsave("UnitRoot_Independent_phi1_bias.eps", device = "eps", width = 6, height = 4, units = "in")

#----Figure 6b, Ly Type 1----
phi1_type1$model <- factor(phi1_type1$model, levels = c("pa", "ds", "gecm", "adl", "adl22", "adl33"), labels = c("pa", "ds", "gecm", "adl", "adl22", "adl33"))
ggplot(phi1_type1, aes(shape = model)) +
  geom_hline(yintercept = 0.05, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = phi1_type1), position = position_dodge(width = 0.65)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 2, byrow = TRUE)) +
  labs(y = "Type 1 Error Rate", x= "T", shape = "", title = expression("Type 1 error rate in the coefficient of y"[t-1])) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("pa", "ds", "gecm", "adl", "adl22", "adl33"),
                     values = c(2, 8, 3, 4, 5, 6),
                     labels = c("PA", 'Dead Start', "GECM", "ADL(1,1)", "ADL(2,2)", "ADL(3,3)")) +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1)
ggsave("UnitRoot_Independent_phi1_type1.pdf", width = 6, height = 4, units = "in")
ggsave("UnitRoot_Independent_phi1_type1.eps", device = "eps", width = 6, height = 4, units = "in")

#----Figure 6c, Ly Power----
phi1_power$model <- factor(phi1_power$model, levels = c("pa", "ds", "gecm", "adl", "adl22", "adl33"), labels = c("pa", "ds", "gecm", "adl", "adl22", "adl33"))
ggplot(phi1_power, aes(shape = model)) +
  geom_hline(yintercept = 1, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = phi1_power), position = position_dodge(width = 0.65)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 2, byrow = TRUE)) +
  labs(y = "Power", x= "T", shape = "", title = expression("Power for the coefficient of y"[t-1])) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("pa", "ds", "gecm", "adl", "adl22", "adl33"),
                     values = c(2, 8, 3, 4, 5, 6),
                     labels = c("PA", 'Dead Start', "GECM", "ADL(1,1)", "ADL(2,2)", "ADL(3,3)")) +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1)
ggsave("UnitRoot_Independent_phi1_power.pdf", width = 6, height = 4, units = "in")
ggsave("UnitRoot_Independent_phi1_power.eps", device = "eps", width = 6, height = 4, units = "in")

#----Figure 7a, LRM Bias----
lrm_bias$lrm_bias[lrm_bias$model == "pa" & lrm_bias$time == 2] <- 19
lrm_bias$model <- factor(lrm_bias$model, levels = c("pa", "gecm", "adl", "adl22", "adl33"), labels = c("pa", "gecm", "adl", "adl22", "adl33"))
ggplot(lrm_bias, aes(shape = model)) +
  geom_hline(yintercept = 0, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = lrm_bias), position = position_dodge(width = 0.75)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 2, byrow = TRUE)) +
  labs(y = "Bias", x= "T", shape = "", title = expression("Bias in the LRM of x"[t])) + ylim(-20, 20) + 
  annotate(geom = "text", x = 2.1, y = 19, label = "1229.44", size = 2) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("pa", "gecm", "adl", "adl22", "adl33"),
                     values = c(2, 3, 4, 5, 6),
                     labels = c("PA", "GECM", "ADL(1,1)", "ADL(2,2)", "ADL(3,3)")) +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1)
ggsave("UnitRoot_Independent_lrm_bias.pdf" , width = 6, height = 4, units = "in")
ggsave("UnitRoot_Independent_lrm_bias.eps", device = "eps", width = 6, height = 4, units = "in")

#----Figure 7b, LRM Type 1----
lrm_type1$model <- factor(lrm_type1$model, levels = c("pa", "gecm", "adl", "adl22", "adl33"), labels = c("pa", "gecm", "adl", "adl22", "adl33"))
ggplot(lrm_type1, aes(shape = model)) +
  geom_hline(yintercept = 0.05, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = lrm_type1), position = position_dodge(width = 0.75)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 2, byrow = TRUE)) +
  labs(y = "Type 1 Error Rate", x= "T", shape = "", title = expression("Type 1 error rate in the LRM of x"[t])) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("pa", "gecm", "adl", "adl22", "adl33"),
                     values = c(2, 3, 4, 5, 6),
                     labels = c("PA", "GECM", "ADL(1,1)", "ADL(2,2)", "ADL(3,3)")) +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1)
ggsave("UnitRoot_Independent_lrm_type1.pdf" , width = 6, height = 4, units = "in")
ggsave("UnitRoot_Independent_lrm_type1.eps", device = "eps", width = 6, height = 4, units = "in")

#---- Appendix Figures Below ----

#----Figure Appendix C-1a, L2y Bias----
phi2_bias$model <- factor(phi2_bias$model, levels = c("adl22", "adl33"), labels = c("adl22", "adl33"))
ggplot(phi2_bias, aes(shape = model)) +
  geom_hline(yintercept = 0, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = phi2_bias), position = position_dodge(width = 0.35)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Bias", x= "T", shape = "", title = expression("Bias in the coefficient of y"[t-2])) + ylim(-0.5, 0.5) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("adl22", "adl33"),
                     values = c(5, 6),
                     labels = c("ADL(2,2)", "ADL(3,3)")) +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1)
ggsave("UnitRoot_Independent_phi2_bias.pdf", width = 6, height = 4, units = "in")

#----Figure Appendix C-1c, L2y Type 1----
phi2_type1$model <- factor(phi2_type1$model, levels = c("adl22", "adl33"), labels = c("adl22", "adl33"))
ggplot(phi2_type1, aes(shape = model)) +
  geom_hline(yintercept = 0.05, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = phi2_type1), position = position_dodge(width = 0.35)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Type 1 Error Rate", x= "T", shape = "", title = expression("Type 1 error rate in the coefficient of y"[t-2])) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("adl22", "adl33"),
                     values = c(5, 6),
                     labels = c("ADL(2,2)", "ADL(3,3)")) +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1)
ggsave("UnitRoot_Independent_phi2_type1.pdf", width = 6, height = 4, units = "in")

#----Figure Appendix C-1b, L3y Bias----
ggplot(phi3_bias, aes(shape = model)) +
  geom_hline(yintercept = 0, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = phi3_bias), position = position_dodge(width = 0.35)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Bias", x= "T", shape = "", title = expression("Bias in the coefficient of y"[t-3])) + ylim(-0.5, 0.5) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10)) +
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("adl33"),
                     values = c(6),
                     labels = c("ADL(3,3)")) 
ggsave("UnitRoot_Independent_phi3_bias.pdf", width = 6, height = 4, units = "in")

#----Figure Appendix C-1d, L3y Type 1----
ggplot(phi3_type1, aes(shape = model)) +
  geom_hline(yintercept = 0.05, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = phi3_type1), position = position_dodge(width = 0.35)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Type 1 Error Rate", x= "T", shape = "", title = expression("Type 1 error rate in the coefficient of y"[t-3])) + ylim(0, 1) + 
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10)) +
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("adl33"),
                     values = c(6),
                     labels = c("ADL(3,3)"))
ggsave("UnitRoot_Independent_phi3_type1.pdf", width = 6, height = 4, units = "in")


#----Figure Appendix C-2a, L2x Bias----
ggplot(b3_bias, aes(shape = model)) +
  geom_hline(yintercept = 0, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = b3_bias), position = position_dodge(width = 0.35)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Bias", x= "T", shape = "", title = expression("Bias in the coefficient of x"[t-2])) + ylim(-0.5, 0.5) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("adl22", "adl33"),
                     values = c(5, 6),
                     labels = c("ADL(2,2)", "ADL(3,3)")) +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1)
ggsave("UnitRoot_Independent_b3_bias.pdf", width = 6, height = 4, units = "in")

#----Figure Appendix C-2c, L2x Type 1----
ggplot(b3_type1, aes(shape = model)) +
  geom_hline(yintercept = 0.05, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = b3_type1), position = position_dodge(width = 0.35)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Type 1 Error Rate", x= "T", shape = "", title = expression("Type 1 error rate in the coefficient of x"[t-2])) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("adl22", "adl33"),
                     values = c(5, 6),
                     labels = c("ADL(2,2)", "ADL(3,3)")) +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1)
ggsave("UnitRoot_Independent_b3_type1.pdf", width = 6, height = 4, units = "in")

#----Figure Appendix C-2b, L3x Bias----
ggplot(b4_bias, aes(shape = model)) +
  geom_hline(yintercept = 0, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = b4_bias), position = position_dodge(width = 0.35)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Bias", x= "T", shape = "", title = expression("Bias in the coefficient of x"[t-3])) + ylim(-0.5, 0.5) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10)) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("adl33"),
                     values = c(6),
                     labels = c("ADL(3,3)")) 
ggsave("UnitRoot_Independent_b4_bias.pdf" , width = 6, height = 4, units = "in")

#----Figure Appendix C-2d, L3x Type 1----
ggplot(b4_type1, aes(shape = model)) +
  geom_hline(yintercept = 0.05, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = b4_type1), position = position_dodge(width = 0.35)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Type 1 Error Rate", x= "T", shape = "", title = expression("Type 1 error rate in the coefficient of x"[t-3])) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10)) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("adl33"),
                     values = c(6),
                     labels = c("ADL(3,3)"))
ggsave("UnitRoot_Independent_b4_type1.pdf" , width = 6, height = 4, units = "in")

#---- Appendix Table Results From Figures Below w/ Additional Coverage Plots ----

#----Figure Appendix C-Table 2, Lags of y Type 1 (sum = 1)----
ftest_y_t1 %>%
  ggplot(aes(shape = model)) +
  geom_hline(yintercept = 0.05, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = ftest_y_t1), position = position_dodge(width = 0.35)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Type 1 Error Rate", x= "T", shape = "", title = expression("Type 1 error rate for a F-test on the coefficients of the lags of y"[t])) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("adl22", "adl33"),
                     values = c(5, 6),
                     labels = c("ADL(2,2)", "ADL(3,3)")) +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1)
ggsave("UnitRoot_Independent_ftest_y_t1.pdf" , width = 6, height = 4, units = "in")

#----Figure Appendix C-Table 2, Lags of y Coverage (sum = 1)----
ftest_y_t1 %>%
  mutate(ftest_y_coverage = 1 - ftest_y_t1) %>%
  ggplot(aes(shape = model)) +
  geom_hline(yintercept = 0.95, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = ftest_y_coverage), position = position_dodge(width = 0.35)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Coverage Probability", x= "T", shape = "", title = expression("Coverage probability for a F-test on the coefficients of the lags of y"[t])) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("adl22", "adl33"),
                     values = c(5, 6),
                     labels = c("ADL(2,2)", "ADL(3,3)")) +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1)
ggsave("UnitRoot_Independent_ftest_y_coverage.pdf" , width = 6, height = 4, units = "in")

#----Figure Appendix C-Table 2, Lags of y Power (sum = 0)----
ftest_y_power %>%
  ggplot(aes(shape = model)) +
  geom_hline(yintercept = 1, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = ftest_y_power), position = position_dodge(width = 0.35)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Power", x= "T", shape = "", title = expression("Power for a F-test on the coefficients of the lags of y"[t])) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("adl22", "adl33"),
                     values = c(5, 6),
                     labels = c("ADL(2,2)", "ADL(3,3)")) +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1)
ggsave("UnitRoot_Independent_ftest_y_power.pdf" , width = 6, height = 4, units = "in")

#----Figure Appendix C-Table 2, Lags of y Power (equal to each other = 0)----
ftest_y_power2 %>%
  ggplot(aes(shape = model)) +
  geom_hline(yintercept = 1, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = ftest_y_power), position = position_dodge(width = 0.35)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Power", x= "T", shape = "", title = expression("Power for a F-test on the coefficients of the lags of y"[t])) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("adl22", "adl33"),
                     values = c(5, 6),
                     labels = c("ADL(2,2)", "ADL(3,3)")) +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1)
ggsave("UnitRoot_Independent_ftest_y_power-2.pdf" , width = 6, height = 4, units = "in")

#----Figure Appendix C-Table 2, x and lags of x Type 1 (sum = 0)----
ftest_x_t1 %>%
  ggplot(aes(shape = model)) +
  geom_hline(yintercept = 0.05, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = ftest_x_t1), position = position_dodge(width = 0.35)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Type 1 Error Rate", x= "T", shape = "", title = expression("Type 1 error rate for a F-test on the coefficients of x"[t]~"and its lags")) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("adl22", "adl33"),
                     values = c(5, 6),
                     labels = c("ADL(2,2)", "ADL(3,3)")) +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1)
ggsave("UnitRoot_Independent_ftest_x_t1.pdf" , width = 6, height = 4, units = "in")

#----Figure Appendix C-Table 2, x and lags of x Coverage (sum = 0)----
ftest_x_t1 %>%
  mutate(ftest_x_coverage = 1 - ftest_x_t1) %>%
  ggplot(aes(shape = model)) +
  geom_hline(yintercept = 0.95, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = ftest_x_coverage), position = position_dodge(width = 0.35)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Coverage Probability", x= "T", shape = "", title = expression("Coverage probability for a F-test on the coefficients of x"[t]~"and its lags")) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("adl22", "adl33"),
                     values = c(5, 6),
                     labels = c("ADL(2,2)", "ADL(3,3)")) +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1)
ggsave("UnitRoot_Independent_ftest_x_coverage.pdf" , width = 6, height = 4, units = "in")

#----Figure Appendix C-Table 2, x and lags of x Type 1 (equal to each other = 0)----
ftest_x_t12 %>%
  ggplot(aes(shape = model)) +
  geom_hline(yintercept = 0.05, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = ftest_x_t1), position = position_dodge(width = 0.35)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Type 1 Error Rate", x= "T", shape = "", title = expression("Type 1 error rate for a F-test on the coefficients of x"[t]~"and its lags")) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("adl22", "adl33"),
                     values = c(5, 6),
                     labels = c("ADL(2,2)", "ADL(3,3)")) +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1)
ggsave("UnitRoot_Independent_ftest_x_t1-2.pdf" , width = 6, height = 4, units = "in")


#----Figure Appendix C-Table 2, x and lags of x Coverage (equal to each other = 0)----
ftest_x_t12 %>%
  mutate(ftest_x_coverage = 1 - ftest_x_t1) %>%
  ggplot(aes(shape = model)) +
  geom_hline(yintercept = 0.95, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = ftest_x_coverage), position = position_dodge(width = 0.35)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Coverage Probability", x= "T", shape = "", title = expression("Coverage probability for a F-test on the coefficients of x"[t]~"and its lags")) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("adl22", "adl33"),
                     values = c(5, 6),
                     labels = c("ADL(2,2)", "ADL(3,3)")) +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1)
ggsave("UnitRoot_Independent_ftest_x_coverage-2.pdf" , width = 6, height = 4, units = "in")


#---- Coverage Plot for Interested Readers Below ----

#----ADL(1,1) model Coverage - x ----
b1_type1 %>% 
  mutate(b1_coverage = 1 - b1_type1) %>%
  filter(model == "adl") %>%
  ggplot() +
  geom_hline(yintercept = 0.95, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = b1_coverage), shape = 4) +
  theme_bw() +
  labs(y = "Coverage Probability", x= "T", title = expression("Coverage probability for the coefficient of x"[t])) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10)) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  theme(legend.position = "none")
ggsave("UnitRoot_Independent_b1_coverage_ADL11.pdf" , width = 6, height = 4, units = "in")

#----ADL(1,1) model Coverage - Lx ----
b2_type1 %>% 
  mutate(b2_coverage = 1 - b2_type1) %>%
  filter(model == "adl") %>%
  ggplot() +
  geom_hline(yintercept = 0.95, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = b2_coverage), shape = 4) +
  theme_bw() +
  labs(y = "Coverage Probability", x= "T", title = expression("Coverage probability for the coefficient of x"[t-1])) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10)) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  theme(legend.position = "none") 
ggsave("UnitRoot_Independent_b2_coverage_ADL11.pdf" , width = 6, height = 4, units = "in")

#----ADL(1,1) model Coverage - Ly ----
phi1_type1 %>%
  mutate(phi1_coverage = 1 - phi1_type1) %>%
  filter(model == "adl") %>%
  ggplot() +
  geom_hline(yintercept = 0.95, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = phi1_coverage), shape = 4) +
  theme_bw() +
  labs(y = "Coverage Probability", x= "T", title = expression("Coverage probability for the coefficient of y"[t-1])) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10)) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  theme(legend.position = "none")
ggsave("UnitRoot_Independent_phi1_coverage_ADL11.pdf" , width = 6, height = 4, units = "in")

#---- Coverage Plot for Corresponding to Discussion in Appendix C2 ----
#----Ly Coverage ----
phi1_type1 %>%
  mutate(phi1_coverage = 1 - phi1_type1) %>%
  ggplot(aes(shape = model)) +
    geom_hline(yintercept = 0.95, color = gray(1/2), lty = 2) +
    geom_point(aes(x = time, y = phi1_coverage), position = position_dodge(width = 0.65)) +
    theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 2, byrow = TRUE)) +
    labs(y = "Coverage Probability", x= "T", shape = "", title = expression("Coverage probability for the coefficient of y"[t-1])) +
    ylim(0, 1) +
    theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
          panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
          panel.grid.minor.y = element_blank()) + 
    scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
    scale_shape_manual(breaks = c("pa", "ds", "gecm", "adl", "adl22", "adl33"),
                       values = c(2, 8, 3, 4, 5, 6),
                       labels = c("PA", 'Dead Start', "GECM", "ADL(1,1)", "ADL(2,2)", "ADL(3,3)")) +
    annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
    annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
    annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
    annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
    annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
    annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
    annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1)
ggsave("UnitRoot_Independent_phi1_coverage.pdf" , width = 6, height = 4, units = "in")

#----L2y Coverage ----
phi2_type1 %>%
  mutate(phi2_coverage = 1 - phi2_type1) %>%
  ggplot(aes(shape = model)) +
    geom_hline(yintercept = 0.95, color = gray(1/2), lty = 2) +
    geom_point(aes(x = time, y = phi2_coverage), position = position_dodge(width = 0.35)) +
    theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
    labs(y = "Coverage Probability", x= "T", shape = "", title = expression("Coverage probability for the coefficient of y"[t-2])) + ylim(0, 1) +
    theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
          panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
          panel.grid.minor.y = element_blank()) + 
    scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
    scale_shape_manual(breaks = c("adl22", "adl33"),
                       values = c(5, 6),
                       labels = c("ADL(2,2)", "ADL(3,3)")) +
    annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
    annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
    annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
    annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
    annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
    annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
    annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1)
ggsave("UnitRoot_Independent_phi2_coverage.pdf" , width = 6, height = 4, units = "in")

#----L3y Coverage ----
phi3_type1 %>%
  mutate(phi3_coverage = 1 - phi3_type1) %>%
ggplot(aes(shape = model)) +
  geom_hline(yintercept = 0.95, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = phi3_coverage), position = position_dodge(width = 0.35)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Coverage Probability", x= "T", shape = "", title = expression("Coverage probability for the coefficient of y"[t-3])) + ylim(0, 1) + 
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10)) +
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("adl33"),
                     values = c(6),
                     labels = c("ADL(3,3)"))
ggsave("UnitRoot_Independent_phi3_coverage.pdf" , width = 6, height = 4, units = "in")

#---- x Coverage ----
b1_type1 %>%
  mutate(b1_coverage = 1 - b1_type1) %>%
  ggplot(aes(shape = model)) +
  geom_hline(yintercept = 0.95, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = b1_coverage), position = position_dodge(width = 0.75)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 2, byrow = TRUE)) +
  labs(y = "Coverage Probability", x= "T", shape = "", title = expression("Coverage probability for the coefficient of x"[t])) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("stat", "diff", "pa", "gecm", "adl", "adl22", "adl33"),
                     values = c(0, 1, 2, 3, 4, 5, 6),
                     labels = c("Static", "Differences", "PA", "GECM", "ADL(1,1)", "ADL(2,2)", "ADL(3,3)")) +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1)
ggsave("UnitRoot_Independent_b1_coverage.pdf" , width = 6, height = 4, units = "in")

#---- Lx Coverage ----
b2_type1 %>%
  mutate(b2_coverage = 1 - b2_type1) %>%
  ggplot(aes(shape = model)) +
  geom_hline(yintercept = 0.95, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = b2_coverage), position = position_dodge(width = 0.75)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 2, byrow = TRUE)) +
  labs(y = "Coverage Probability", x= "T", shape = "", title = expression("Coverage probability for the coefficient of x"[t-1])) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("ds", "gecm", "adl", "adl22", "adl33"),
                     values = c(8, 3, 4, 5, 6),
                     labels = c("Dead Start", "GECM", "ADL(1,1)", "ADL(2,2)", "ADL(3,3)")) +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1)
ggsave("UnitRoot_Independent_b2_coverage.pdf" , width = 6, height = 4, units = "in")

#---- L2x Coverage ----
b3_type1 %>%
  mutate(b3_coverage = 1 - b3_type1) %>%
  ggplot(aes(shape = model)) +
  geom_hline(yintercept = 0.95, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = b3_coverage), position = position_dodge(width = 0.35)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Coverage Probability", x= "T", shape = "", title = expression("Coverage probability the coefficient of x"[t-2])) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("adl22", "adl33"),
                     values = c(5, 6),
                     labels = c("ADL(2,2)", "ADL(3,3)")) +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1)
ggsave("UnitRoot_Independent_b3_coverage.pdf" , width = 6, height = 4, units = "in")

#---- L3x Coverage ----
b4_type1 %>%
  mutate(b4_coverage = 1 - b4_type1) %>%
  ggplot(aes(shape = model)) +
  geom_hline(yintercept = 0.95, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = b4_coverage), position = position_dodge(width = 0.35)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 1, byrow = TRUE)) +
  labs(y = "Coverage Probability", x= "T", shape = "", title = expression("Coverage probability for the coefficient of x"[t-3])) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10)) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("adl33"),
                     values = c(6),
                     labels = c("ADL(3,3)"))
ggsave("UnitRoot_Independent_b4_coverage.pdf" , width = 6, height = 4, units = "in")

#---- LRM Coverage ----
lrm_type1 %>%
  mutate(lrm_coverage = 1 - lrm_type1) %>%
ggplot(aes(shape = model)) +
  geom_hline(yintercept = 0.95, color = gray(1/2), lty = 2) +
  geom_point(aes(x = time, y = lrm_coverage), position = position_dodge(width = 0.75)) +
  theme_bw() + theme(legend.position = "bottom") + guides(shape = guide_legend(nrow = 2, byrow = TRUE)) +
  labs(y = "Coverage Probability", x= "T", shape = "", title = expression("Coverage probability the LRM of x"[t])) + ylim(0, 1) +
  theme(plot.title = element_text(size = 10), axis.title = element_text(size = 10),
        panel.grid.major.x = element_blank(), axis.ticks.x = element_blank(),
        panel.grid.minor.y = element_blank()) + 
  scale_x_discrete(limits = c("25", "50", "75", "100", "150", "500", "1000")) +
  scale_shape_manual(breaks = c("pa", "gecm", "adl", "adl22", "adl33"),
                     values = c(2, 3, 4, 5, 6),
                     labels = c("PA", "GECM", "ADL(1,1)", "ADL(2,2)", "ADL(3,3)")) +
  annotate("rect", xmin = 0.65, xmax = 1.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 1.65, xmax = 2.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 2.65, xmax = 3.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 3.65, xmax = 4.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 4.65, xmax = 5.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 5.65, xmax = 6.35, ymin = -Inf, ymax = Inf, alpha = .1) +
  annotate("rect", xmin = 6.65, xmax = 7.35, ymin = -Inf, ymax = Inf, alpha = .1)
ggsave("UnitRoot_Independent_lrm_coverage.pdf" , width = 6, height = 4, units = "in")

